Detailed balance #52

henrikjacobsenfys · 2025-09-23T09:57:49Z

Utility function to calculate the detailed balance factor as discussed here: easyscience/dynamics#6

Includes the option to normalise by temperature. Takes various inputs and outputs a np array. I may want to later output a scipp array, but not yet.

codecov · 2025-09-23T10:01:42Z

Codecov Report

❌ Patch coverage is 91.35802% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.35%. Comparing base (3d99faa) to head (550bbc6).
⚠️ Report is 1 commits behind head on develop.

Files with missing lines	Patch %	Lines
src/easydynamics/utils/detailed_balance.py	91.13%	7 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff              @@
##           develop      #52       +/-   ##
============================================
+ Coverage     0.00%   91.35%   +91.35%     
============================================
  Files            1        2        +1     
  Lines            2       81       +79     
============================================
+ Hits             0       74       +74     
- Misses           2        7        +5

Flag	Coverage Δ
unittests	`91.35% <91.35%> (+91.35%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Co-authored-by: rozyczko <[email protected]>

rozyczko

Good idea with separating utilities from the rest of the code.
A few points added on how to improve the code.
The rest looks good!

src/easydynamics/utils/utils.py

rozyczko · 2025-09-23T11:50:59Z

src/easydynamics/utils/utils.py

+    """
+
+    # Input validation
+    if not isinstance(divide_by_temperature, bool):


Type checks already test the types. If we enforce those, these checks aren't necessary

If this is a user-facing function, input validation is always a good idea since you have no guarantee that the users have any clue what they're doing :)

So it's mostly to be used internally, but I can't rule out that users will want to use it.

You should be perfectly clear wether a function/method is user-facing or not.
For user-facing functions you need input validation checks like these and thoroughly documented docstrings.
For developer/hidden functions you don't as you always know what you pass to the function. Also, those methods should always start with an _.

src/easydynamics/utils/utils.py

henrikjacobsenfys · 2025-09-24T13:47:41Z

@rozyczko Just double checking, is it OK to merge? :)

rozyczko

looks good now!

damskii9992

Will continue tomorrow

damskii9992 · 2025-10-02T13:43:35Z

src/easydynamics/utils/utils.py

+    """
+
+    # Input validation
+    if not isinstance(divide_by_temperature, bool):


If this is a user-facing function, input validation is always a good idea since you have no guarantee that the users have any clue what they're doing :)

damskii9992 · 2025-10-02T13:44:23Z

src/easydynamics/utils/utils.py

+
+    # Input validation
+    if not isinstance(divide_by_temperature, bool):
+        raise ValueError("divide_by_temperature must be a boolean.")


I would write "must be True or False" rather than "must be boolean", slightly easier to understand for non-programmers ;)

Also, this should be a type error

damskii9992 · 2025-10-02T13:46:59Z

src/easydynamics/utils/utils.py

+    # Convert to numpy array first for consistent handling
+    if isinstance(value, (int, float)):
+        array_value = np.array(value)
+    elif isinstance(value, (list, tuple)):


Do you also allow tuples? That is not shown in the type hints

damskii9992 · 2025-10-02T13:51:21Z

src/easydynamics/utils/utils.py

+        DBF_values = DBF.values
+    return DBF_values
+
+def _convert_to_scipp_variable(value: Union[int, float, list, np.ndarray, Parameter, sc.Variable], unit: str, name: str = "value") -> sc.Variable:


Why is unit a mandatory field? Especially since you reuse the unit if the input is a parameter?

damskii9992 · 2025-10-02T13:52:23Z

src/easydynamics/utils/utils.py

+        return sc.scalar(value=float(array_value.flat[0]), unit=unit)
+    else:
+        # Multi-element array
+        return sc.array(dims=['x'], values=array_value, unit=unit)


Is this method meant to take multi-dimensional lists too? Won't you need more than 1 dim then?

Not at the moment, no. It may also be that I need to revisit the dims if/when I move entirely to scipp for x, but I'll park that with the other headaches from the other PR for now :)

damskii9992 · 2025-10-02T13:55:22Z

src/easydynamics/utils/utils.py

+        raise ValueError("energy_unit must be a string.")
+
+    if not isinstance(temperature_unit, str):
+        raise ValueError("temperature_unit must be a string.")


damskii9992 · 2025-10-02T14:00:55Z

src/easydynamics/utils/utils.py

+        raise ValueError("Temperature must be non-negative.")
+
+    # Convert energy to sc variable to make units easy to handle
+    energy = _convert_to_scipp_variable(energy, energy_unit, "energy")


The error message here would say that energy can be a Parameter, which it can't according to your type hints.

damskii9992 · 2025-10-02T14:06:43Z

src/easydynamics/utils/utils.py

+    # We give users the option to specify the unit of the energy, but if the input has a unit, they might clash
+    if energy.unit != energy_unit:
+        warnings.warn(f"Input energy has unit {energy.unit}, but energy_unit was set to {energy_unit}. Using {energy.unit}.")
+        energy_unit = energy.unit


I see what you're doing here. I think I would do this slightly differently. Rather than having an energy_unit at all, I would simply raise a warning stating that the input was assumed to be in meV if it has no units. If the user just blindly passes a numpy array right now, they would not be informed that it was assumed to be in meV automatically.

They're informed in the docstring - is that enough?
Args:
energy : number, list, np.ndarray, or scipp Variable. If number, assumed to be in meV unless energy_unit is set.
Energy transfer

damskii9992 · 2025-10-02T14:07:15Z

src/easydynamics/utils/utils.py

+    # Same for temperature
+    if temperature.unit != temperature_unit:
+        warnings.warn(f"Input temperature has unit {temperature.unit}, but temperature_unit was set to {temperature_unit}. Using {temperature.unit}.")
+        temperature_unit = temperature.unit


And again, the same here.

damskii9992 · 2025-10-02T14:08:18Z

src/easydynamics/utils/utils.py

+    # Zero temperature deserves special treatment
+    if temperature.value == 0:
+        if divide_by_temperature:
+            raise ValueError("Cannot divide by T when T=0.")


You can raise a specific ZeroDivisionError here ;)

damskii9992

More comments

damskii9992 · 2025-10-03T09:14:36Z

src/easydynamics/utils/detailed_balance.py

+
+def _convert_to_scipp_variable(
+    value: Union[int, float, list, np.ndarray, Parameter, sc.Variable],
+    unit: Optional[str],


Type hinting the unit argument does not make it optional :)
You need to add a default value for it, like None, if you wish to make it optional.

damskii9992 · 2025-10-03T09:56:54Z

src/easydynamics/utils/detailed_balance.py

+        # Scalar or single-element array
+        try:
+            return sc.scalar(value=float(array_value.flat[0]), unit=unit)
+        except sc.UnitError as e:


Since you imported the UnitError from scipp, you don't need the "sc." prefix here or below :)

damskii9992 · 2025-10-03T12:32:11Z

src/easydynamics/utils/detailed_balance.py

+            energy < 0.0 * sc.Unit(energy_unit), 0.0 * sc.Unit(energy_unit), energy
+        )
+
+        if DBF.sizes == {}:


. . . . Does scipp return an empty dictionary when you get sizes on an empty array? O.o

damskii9992 · 2025-10-03T12:38:04Z

src/easydynamics/utils/detailed_balance.py

+    # Small x (small energy and/or high temperature): Taylor expansion
+    small = sc.abs(x) < SMALL_THRESHOLD
+
+    DBF = sc.where(small, 1 + x / 2 + x**2 / 12, DBF)


List comprehension

damskii9992 · 2025-10-03T12:46:54Z

src/easydynamics/utils/detailed_balance.py

+    DBF = sc.where(large, x, DBF)
+
+    # General case: exact formula
+    mid = ~small & ~large


~ should be not

nvm. remove tildes and replace & with ==

damskii9992 · 2025-10-03T12:48:39Z

src/easydynamics/utils/detailed_balance.py

+
+    # General case: exact formula
+    mid = ~small & ~large
+    DBF = sc.where(mid, x / (1 - sc.exp(-x)), DBF)


Write comment: zeros handles by small treshold

seventil · 2025-10-03T13:18:33Z

src/easydynamics/utils/detailed_balance.py

+            Temperature
+        energy_unit : str, optional
+            Unit for energy if energy is given as a number or list. Default is 'meV'
+        temperature_unit : str, optional


it refers to the formula instead of the function parameter

seventil · 2025-10-03T13:26:56Z

src/easydynamics/utils/detailed_balance.py

+        warnings.warn(
+            f"Input energy has unit {energy.unit}, but energy_unit was set to {energy_unit}. Using {energy.unit}."
+        )
+        energy_unit = energy.unit


why not use energy.unit in code below after converting to sc? That way, energy_unit will only be an argument for convert function and there won't be two variables with potentially different contents referring to essentially the same entity. Same for temperature

seventil · 2025-10-03T13:31:34Z

src/easydynamics/utils/detailed_balance.py

+    # Large x (large energy and/or low temperature): asymptotic form
+    large = x > LARGE_THRESHOLD
+    DBF = sc.where(large, x, DBF)
+


why is x in small casted to abs, but not in large?

seventil · 2025-10-03T13:35:39Z

src/easydynamics/utils/detailed_balance.py

+        array_value = np.array(value.value)
+        unit = value.unit
+    else:
+        if name == "energy":


if value is Parameter and name is "energy" it would be checked out by the elif isinstance(value, Parameter): in the upper part, this specific if statement for "energy" avoiding Parameter type does not make sense

seventil · 2025-10-03T13:44:35Z

src/easydynamics/utils/detailed_balance.py

+def _convert_to_scipp_variable(
+    value: Union[int, float, list, np.ndarray, Parameter, sc.Variable],
+    unit: Optional[str] = None,
+    name: str = "value",


the default value is not used anywhere except error messages and it's not marked as optional.
imho it's a redundant parameter since you'll see the variable name that calls _convert_to_scipp_variable in the traceback anyway

damskii9992

Final comments to the tests now. The code itself looks good.

damskii9992 · 2025-10-06T09:31:59Z

tests/unit_tests/utils/test_detailed_balance.py

+from easyscience import Parameter
+from scipp.constants import Boltzmann as kB
+
+from easydynamics.utils import _detailed_balance_factor as detailed_balance_factor


Don't remove the underscore in imports. At least not in code.

damskii9992 · 2025-10-06T09:35:27Z

tests/unit_tests/utils/test_detailed_balance.py

+        temperature_unit = 5
+        # Then Expect
+        with pytest.raises(TypeError, match="temperature_unit must be a string."):
+            detailed_balance_factor(energy, T, temperature_unit=temperature_unit)


It is a good idea to test multiple wrong inputs, not just a single one (here, you chose an int)

damskii9992 · 2025-10-06T09:58:22Z

tests/unit_tests/utils/test_detailed_balance.py

+    def test_scalar_input(self):
+        # When
+        energy = 2.0
+        T = 100


Considering how often you set energy and temperature in your # When, it might be worth making this as a fixture to reuse throughout your tests.

damskii9992 · 2025-10-06T10:01:55Z

tests/unit_tests/utils/test_detailed_balance.py

+            energy / (1 - np.exp(-energy / (kB_meV_per_K * T))) / (kB_meV_per_K * T)
+        )
+
+        np.testing.assert_allclose(result, expected, rtol=1e-5)


Since you're testing for a scalar input, it probably makes sense to test that the output is also a scalar? :)

It's an array, but good point

damskii9992 · 2025-10-06T10:06:46Z

tests/unit_tests/utils/test_detailed_balance.py

+
+        assert isinstance(result, np.ndarray)
+        assert result.shape == (3,)
+        np.testing.assert_allclose(result, expected_values, rtol=1e-5)


All of these tests could have been a single parametrized test with the energy input parametrized

I will keep this in mind for future tests ;)

damskii9992 · 2025-10-06T10:10:19Z

tests/unit_tests/utils/test_detailed_balance.py

+    def test_negative_temperature_raises(self):
+        # When Then Expect
+        with pytest.raises(ValueError, match="Temperature must be non-negative"):
+            detailed_balance_factor(1.0, -10)


The temperature tests can be collected and parametrized and the exception tests can also be collected and parametrized.

damskii9992 · 2025-10-06T12:10:06Z

tests/unit_tests/utils/test_utils.py

+        )
+        # Expect
+        expected = np.full(5, kB_meV_per_K * T)
+        np.testing.assert_allclose(result, expected, rtol=1e-5)


Shouldn't the expected value here also change with energy? Why are you simply comparing to a constant?
The way you're testing it right now, with a relative tolerance equal to your highest energy, you're not really testing anything.

Good find, I messed that test up.

damskii9992 · 2025-10-06T12:12:01Z

tests/unit_tests/utils/test_utils.py

+            energy=energy, temperature=T, divide_by_temperature=False
+        )
+        # Expect
+        np.testing.assert_allclose(result, energy, rtol=1e-2)


I'm also not sure this test makes much sense. In the high-energy limit you set DBF = x. So shouldn't you just check that they're literally equal here? Instead of testing that they're close?

They differ on like 13th decimal, probably due to floating point things. I made the test much stricter.

damskii9992 · 2025-10-06T12:19:23Z

tests/unit_tests/utils/test_utils.py

+
+        # Expect
+        expected_ratio = np.exp(energy / (kB_meV_per_K * T))
+        np.testing.assert_allclose(ratio, expected_ratio, rtol=1e-5)


Maybe add a comment for the reasoning behind this test? It is not obvious why the ratios should fulfill this condition. Think about your future programmers :)

But everyone knows the principle of detailed balance, right? (I've added a comment)

damskii9992 · 2025-10-06T12:25:18Z

tests/unit_tests/utils/test_utils.py

+        expected = energy.values / (
+            1 - np.exp(-energy.values / 1000 / (kB_meV_per_K * T))
+        )
+        np.testing.assert_allclose(result, expected, rtol=1e-5)


You don't have a test of a faulty unit such as "m" or "l"

In both this and the temperature unit. I also don't think you have a test where you pass a sc.Unit?

henrikjacobsenfys added 2 commits September 23, 2025 11:45

Add detailed balance and tests

7767fe3

add __init__

64453d5

Add missing test

7bf99c9

Copilot AI mentioned this pull request Sep 23, 2025

Add comprehensive critique of PR #52 "Detailed balance" #53

Closed

Copilot AI added a commit that referenced this pull request Sep 23, 2025

Complete comprehensive critique of PR #52 detailed balance

29e2713

Co-authored-by: rozyczko <[email protected]>

Add another test

cb47df9

henrikjacobsenfys added [scope] enhancement Adds/improves features (major.MINOR.patch) [priority] medium Normal/default priority labels Sep 23, 2025

Make copilot happy

87353c0

rozyczko reviewed Sep 23, 2025

View reviewed changes

henrikjacobsenfys added 2 commits September 23, 2025 15:08

Fixing minor issues

b9831ba

Update example and thresholds

c036a97

rozyczko approved these changes Oct 2, 2025

View reviewed changes

damskii9992 requested changes Oct 2, 2025

View reviewed changes

Add some checks

18f88dd

damskii9992 reviewed Oct 3, 2025

View reviewed changes

Fix various small things

3bbaed3

seventil reviewed Oct 3, 2025

View reviewed changes

henrikjacobsenfys added 3 commits October 4, 2025 20:24

Minor improvements, updated .gitignore

349fe04

Fix unit error

ba1d543

Add test of incompatible units

99f2e12

damskii9992 requested changes Oct 6, 2025

View reviewed changes

henrikjacobsenfys added 2 commits October 6, 2025 15:37

Final touches

760c729

final touches

550bbc6

damskii9992 approved these changes Oct 7, 2025

View reviewed changes

henrikjacobsenfys merged commit d51e07e into develop Oct 7, 2025
16 checks passed

Detailed balance #52

Detailed balance #52

Uh oh!

Conversation

henrikjacobsenfys commented Sep 23, 2025

Uh oh!

codecov bot commented Sep 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

rozyczko left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

henrikjacobsenfys commented Sep 24, 2025

Uh oh!

rozyczko left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

damskii9992 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

damskii9992 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Sep 23, 2025 •

edited

Loading

rozyczko left a comment •

edited

Loading